! This file is part of xtb.
!
! Copyright (C) 2019-2020 Stefan Grimme
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb.  If not, see <https://www.gnu.org/licenses/>.

      subroutine gfnffdrab(n,at,cn,dcn,nsrb,srblist,rab,grab,rabdcn)
      use xtb_mctc_accuracy, only : wp
      implicit none
      integer n                 ! number of atoms
      integer at(n)             ! ordinal numbers
      real*8  cn(n)             ! D3 CN
      real*8  dcn(3,n,n)        ! D3 CN derivatives
      integer nsrb              ! # of pairs
      integer srblist(3,nsrb)   ! list of atom pairs
      real*8 rab(nsrb)          ! output bond lengths estimates, input the predetermined bond shifts
      real*8 grab(3,n,nsrb)     ! output bond lengths gradients
      real*8 rabdcn(2,nsrb)       ! output bond lengths derivative w.r.t. cn

!     local variables
      integer m,i,j,k,ii,jj,ati,atj,ir,jr
      interface
         integer function itabrow6(i)
            integer i
         end function
      end interface

      real*8 cnfak(103),r0(103),en(103)
      real*8 ra,rb,k1,k2,den,ff
      real(wp) :: p(6,2)
      real*8 scaleF(103,103)

c     fitted on PBExa-3c geom set by SG, 9/2018
c     H B C N O F SI P S CL GE AS SE BR SN SB TE I together (and for glob par)
c     and rest one element at a time, LN taken as av. of La and Hf
c     about 15-20 reference molecules per element (about 30-40 data points per fit
c     parameter)
! electronegativity (started from Pauling values)
! base radius
! CN dep. lin
! CN dep. quadratic (not used, only slightly better)

      data en /
     . 2.30085633, 2.78445145, 1.52956084, 1.51714704, 2.20568300,
     . 2.49640820, 2.81007174, 4.51078438, 4.67476223, 3.29383610,
     . 2.84505365, 2.20047950, 2.31739628, 2.03636974, 1.97558064,
     . 2.13446570, 2.91638164, 1.54098156, 2.91656301, 2.26312147,
     . 2.25621439, 1.32628677, 2.27050569, 1.86790977, 2.44759456,
     . 2.49480042, 2.91545568, 3.25897750, 2.68723778, 1.86132251,
     . 2.01200832, 1.97030722, 1.95495427, 2.68920990, 2.84503857,
     . 2.61591858, 2.64188286, 2.28442252, 1.33011187, 1.19809388,
     . 1.89181390, 2.40186898, 1.89282464, 3.09963488, 2.50677823,
     . 2.61196704, 2.09943450, 2.66930105, 1.78349472, 2.09634533,
     . 2.00028974, 1.99869908, 2.59072029, 2.54497829, 2.52387890,
     . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000, ! 60
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, ! 70
     . 2.00000000, 2.30089349, 1.75039077, 1.51785130, 2.62972945,
     . 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898, ! 80
     . 1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095,
     . 2.34224386,
     . 2.52387890,
     . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000, ! 92
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, ! 102
     . 2.00000000
     ./
      data r0 /
     . 0.55682207, 0.80966997, 2.49092101, 1.91705642, 1.35974851,
     . 0.98310699, 0.98423007, 0.76716063, 1.06139799, 1.17736822, ! 10
     . 2.85570926, 2.56149012, 2.31673425, 2.03181740, 1.82568535,
     . 1.73685958, 1.97498207, 2.00136196, 3.58772537, 2.68096221, ! 20
     . 2.23355957, 2.33135502, 2.15870365, 2.10522128, 2.16376162,
     . 2.10804037, 1.96460045, 2.00476257, 2.22628712, 2.43846700, ! 30
     . 2.39408483, 2.24245792, 2.05751204, 2.15427677, 2.27191920,
     . 2.19722638, 3.80910350, 3.26020971, 2.99716916, 2.71707818, ! 40
     . 2.34950167, 2.11644818, 2.47180659, 2.32198800, 2.32809515,
     . 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278, ! 50
     . 2.56912355, 2.54374096, 2.56914830, 2.53680807, 4.24537037,
     . 3.66542289, 3.39000000, 3.22000000, 3.19000000, 3.17000000, ! 60
     . 3.21000000, 3.05000000, 3.08000000, 3.08000000, 3.06000000,
     . 2.92000000, 2.99000000, 3.00000000, 2.91000000, 2.92000000, ! 70
     . 2.76000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577,
     . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, ! 80
     . 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552,
     . 2.63264944,
     . 4.24537037,
     . 3.66542289,3.440000000,3.110000000,2.750000000,2.690000000,
     . 2.97000000,2.860000000,2.860000000,3.010000000,3.180000000,
     . 3.12000000,2.620000000,2.670000000,2.670000000,2.750000000,
     . 2.62000000
     ./
      data cnfak /
     . 0.17957827, 0.25584045,-0.02485871, 0.00374217, 0.05646607,
     . 0.10514203, 0.09753494, 0.30470380, 0.23261783, 0.36752208,
     . 0.00131819,-0.00368122,-0.01364510, 0.04265789, 0.07583916,
     . 0.08973207,-0.00589677, 0.13689929,-0.01861307, 0.11061699,
     . 0.10201137, 0.05426229, 0.06014681, 0.05667719, 0.02992924,
     . 0.03764312, 0.06140790, 0.08563465, 0.03707679, 0.03053526,
     .-0.00843454, 0.01887497, 0.06876354, 0.01370795,-0.01129196,
     . 0.07226529, 0.01005367, 0.01541506, 0.05301365, 0.07066571,
     . 0.07637611, 0.07873977, 0.02997732, 0.04745400, 0.04582912,
     . 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441, ! 50
     . 0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969,
     .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 60
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 70
     . 0.07700000, 0.08379100, 0.07314553, 0.05318438, 0.06799334,
     . 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502, ! 80
     . 0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675,
     . 0.04800662, ! 86
     .-0.04110969, ! 87
     .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 92
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 102
     . 0.07700000 ! 103
     ./

      ! setup factor for Fluorine-Ln or F-Ac bonds
      scaleF = 1.0
      scaleF(57:71,9) =0.93533568
      scaleF(89:103,9)=0.93533568
      scaleF(9,57:71) =0.93533568
      scaleF(9,89:103)=0.93533568
      ! for Chlorine
      scaleF(57:71,17) =1.0190114
      scaleF(89:103,17)=1.0190114
      scaleF(17,57:71) =1.0190114
      scaleF(17,89:103)=1.0190114
      ! for Bromine
      scaleF(57:71,35) =1.0425532
      scaleF(89:103,35)=1.0425532
      scaleF(35,57:71) =1.0425532
      scaleF(35,89:103)=1.0425532
      ! for Iodine
      scaleF(57:71,53) =1.0551948
      scaleF(89:103,53)=1.0551948
      scaleF(53,57:71) =1.0551948
      scaleF(53,89:103)=1.0551948

c     global EN polynominal parameters x 10^3
      p(1,1)=    29.84522887
      p(2,1)=    -1.70549806
      p(3,1)=     6.54013762
      p(4,1)=     6.39169003
      p(5,1)=     6.00000000 ! handish downward extrapolated for GFN-FF
      p(6,1)=     5.60000000
      p(1,2)=    -8.87843763
      p(2,2)=     2.10878369
      p(3,2)=     0.08009374
      p(4,2)=    -0.85808076
      p(5,2)=    -1.15000000
      p(6,2)=    -1.30000000


      do k=1,nsrb
         j = srblist(1,k)
         i = srblist(2,k)
         ati=at(i)
         atj=at(j)
         ir=itabrow6(ati)
         jr=itabrow6(atj)
c--------
         ra=r0(ati)+cnfak(ati)*cn(i)
         rb=r0(atj)+cnfak(atj)*cn(j)
         den=abs(en(ati)-en(atj))
         k1=0.005d0*(p(ir,1)+p(jr,1))
         k2=0.005d0*(p(ir,2)+p(jr,2))
         ff=1.0d0-k1*den-k2*den**2
         rab(k) =(ra+rb+rab(k))*ff*scaleF(ati,atj)
         rabdcn(1,k)=cnfak(ati)*ff*scaleF(ati,atj)
         rabdcn(2,k)=cnfak(atj)*ff*scaleF(ati,atj)
         do m=1,n
            grab(1:3,m,k)=scaleF(ati,atj)*ff*(cnfak(ati)*dcn(1:3,m,i)
     .                       +cnfak(atj)*dcn(1:3,m,j))
         enddo
c--------
      enddo

      end subroutine gfnffdrab

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      INTEGER FUNCTION iTabRow6(i)
      implicit none
      INTEGER i

      iTabRow6=0
      If (i.gt. 0 .and. i.le. 2) Then
         iTabRow6=1
      Else If (i.gt. 2 .and. i.le.10) Then
         iTabRow6=2
      Else If (i.gt.10 .and. i.le.18) Then
         iTabRow6=3
      Else If (i.gt.18 .and. i.le.36) Then
         iTabRow6=4
      Else If (i.gt.36 .and. i.le.54) Then
         iTabRow6=5
      Else If (i.gt.54) Then
         iTabRow6=6
      End If

      Return
      End FUNCTION iTabRow6

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine gfnffrab(n,at,cn,rab)
      implicit none
      integer n                 ! number of atoms
      integer at(n)             ! ordinal numbers
      real*8  cn(n)             ! D3 CN
      real*8 rab(n*(n+1)/2)     ! output bond lengths estimates

      integer m,i,j,k,ii,jj,ati,atj,ir,jr
      INTEGER iTabRow6,lin

      real*8 cnfak(103),r0(103),en(103)
      real*8 ra,rb,k1,k2,den,ff,p(6,2)
      real*8 scaleF(103,103)
      data en /
     . 2.30085633, 2.78445145, 1.52956084, 1.51714704, 2.20568300,
     . 2.49640820, 2.81007174, 4.51078438, 4.67476223, 3.29383610,
     . 2.84505365, 2.20047950, 2.31739628, 2.03636974, 1.97558064,
     . 2.13446570, 2.91638164, 1.54098156, 2.91656301, 2.26312147,
     . 2.25621439, 1.32628677, 2.27050569, 1.86790977, 2.44759456,
     . 2.49480042, 2.91545568, 3.25897750, 2.68723778, 1.86132251,
     . 2.01200832, 1.97030722, 1.95495427, 2.68920990, 2.84503857,
     . 2.61591858, 2.64188286, 2.28442252, 1.33011187, 1.19809388,
     . 1.89181390, 2.40186898, 1.89282464, 3.09963488, 2.50677823,
     . 2.61196704, 2.09943450, 2.66930105, 1.78349472, 2.09634533,
     . 2.00028974, 1.99869908, 2.59072029, 2.54497829, 2.52387890,
     . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.30089349, 1.75039077, 1.51785130, 2.62972945,
     . 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898,
     . 1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095,
     . 2.34224386,
     . 2.52387890,
     . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000, ! 92
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,
     . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, ! 102
     . 2.00000000
     ./
      data r0 /
     . 0.55682207, 0.80966997, 2.49092101, 1.91705642, 1.35974851,
     . 0.98310699, 0.98423007, 0.76716063, 1.06139799, 1.17736822, ! 10
     . 2.85570926, 2.56149012, 2.31673425, 2.03181740, 1.82568535,
     . 1.73685958, 1.97498207, 2.00136196, 3.58772537, 2.68096221, ! 20
     . 2.23355957, 2.33135502, 2.15870365, 2.10522128, 2.16376162,
     . 2.10804037, 1.96460045, 2.00476257, 2.22628712, 2.43846700, ! 30
     . 2.39408483, 2.24245792, 2.05751204, 2.15427677, 2.27191920,
     . 2.19722638, 3.80910350, 3.26020971, 2.99716916, 2.71707818, ! 40
     . 2.34950167, 2.11644818, 2.47180659, 2.32198800, 2.32809515,
     . 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278, ! 50
     . 2.56912355, 2.54374096, 2.56914830, 2.53680807, 4.24537037,
     . 3.66542289, 3.39000000, 3.22000000, 3.19000000, 3.17000000, ! 60
     . 3.21000000, 3.05000000, 3.08000000, 3.08000000, 3.06000000,
     . 2.92000000, 2.99000000, 3.00000000, 2.91000000, 2.92000000, ! 70
     . 2.76000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577,
     . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, ! 80
     . 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552,
     . 2.63264944,
     . 4.24537037,
     . 3.66542289,3.440000000,3.110000000,2.750000000,2.690000000,
     . 2.97000000,2.860000000,2.860000000,3.010000000,3.180000000,
     . 3.12000000,2.620000000,2.670000000,2.670000000,2.750000000,
     . 2.62000000
     ./
      data cnfak /
     . 0.17957827, 0.25584045,-0.02485871, 0.00374217, 0.05646607,
     . 0.10514203, 0.09753494, 0.30470380, 0.23261783, 0.36752208,
     . 0.00131819,-0.00368122,-0.01364510, 0.04265789, 0.07583916,
     . 0.08973207,-0.00589677, 0.13689929,-0.01861307, 0.11061699,
     . 0.10201137, 0.05426229, 0.06014681, 0.05667719, 0.02992924,
     . 0.03764312, 0.06140790, 0.08563465, 0.03707679, 0.03053526,
     .-0.00843454, 0.01887497, 0.06876354, 0.01370795,-0.01129196,
     . 0.07226529, 0.01005367, 0.01541506, 0.05301365, 0.07066571,
     . 0.07637611, 0.07873977, 0.02997732, 0.04745400, 0.04582912,
     . 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441,
     . 0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969,
     .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.08379100, 0.07314553, 0.05318438, 0.06799334,
     . 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502,
     . 0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675,
     . 0.04800662,
     .-0.04110969, ! 87
     .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 92
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,
     . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 102
     . 0.07700000 ! 103
     ./

      ! setup factor for Fluorine-Ln or F-Ac bonds
      scaleF = 1.0
      scaleF(57:71,9) =0.93533568
      scaleF(89:103,9)=0.93533568
      scaleF(9,57:71) =0.93533568
      scaleF(9,89:103)=0.93533568
      ! for Chlorine
      scaleF(57:71,17) =1.0190114
      scaleF(89:103,17)=1.0190114
      scaleF(17,57:71) =1.0190114
      scaleF(17,89:103)=1.0190114
      ! for Bromine
      scaleF(57:71,35) =1.0425532
      scaleF(89:103,35)=1.0425532
      scaleF(35,57:71) =1.0425532
      scaleF(35,89:103)=1.0425532
      ! for Iodine
      scaleF(57:71,53) =1.0551948
      scaleF(89:103,53)=1.0551948
      scaleF(53,57:71) =1.0551948
      scaleF(53,89:103)=1.0551948

      p(1,1)=    29.84522887
      p(2,1)=    -1.70549806
      p(3,1)=     6.54013762
      p(4,1)=     6.39169003
      p(5,1)=     6.00000000
      p(6,1)=     5.60000000
      p(1,2)=    -8.87843763
      p(2,2)=     2.10878369
      p(3,2)=     0.08009374
      p(4,2)=    -0.85808076
      p(5,2)=    -1.15000000
      p(6,2)=    -1.30000000

      do i=1,n
         do j=1,i
         k=lin(j,i)
         ati=at(i)
         atj=at(j)
         ir=itabrow6(ati)
         jr=itabrow6(atj)
         ra=r0(ati)+cnfak(ati)*cn(i)
         rb=r0(atj)+cnfak(atj)*cn(j)
         den=abs(en(ati)-en(atj))
         k1=0.005d0*(p(ir,1)+p(jr,1))
         k2=0.005d0*(p(ir,2)+p(jr,2))
         ff=1.0d0-k1*den-k2*den**2
         rab(k)=(ra+rb)*ff*scaleF(ati,atj)
         enddo
      enddo

      end subroutine gfnffrab
